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ABSTRACT 


A finite element method was applied to the unsteady 
transonic small disturbance equation and integrated until 
the solution converged to the steady state for a thin non- 
lifting airfoil. The method of weighted residuals was used 
to formulate the finite element equations, and Houbolt's 
method of central differencing in time was used to integrate 
these equations. 


A secondary investigation applied the steady transonic 


small disturbance equations to a converging-diverging nozzle. 
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I. INTRODUCTION 


Transonic inviscid flows past a smooth airfoil may be 
expressed in terms of the velocity potential ® satisfying 


the transonic small disturbance equation, 


(LM M20 * LODO + by 70 (1) 


This equation presents two major difficulties, 1) it is non- 
linear and 2) it is of the mixed hyperbolic-elliptic type. 
Analytical solutions to non-linear equations are difficult 
to obtain. One must normally resort to numerical methods. 
When the coefficient (1 - M2 - M2 ty “* 1)o,) in equation 1 
is negative, the flow is supersonic and the equation is 
called hyperbolic; otherwise the flow is subsonic and the 
equation is elliptic. The forms of the two solutions are 
fundamentally different. Hyperbolic equations admit both 
discontinuities, which propagate only in characteristic 
directions, and the presence of shock fronts. Elliptic 
equations, on the other hand, require that the dependent 
variables and their derivatives be continuous and that a 
change in any part of the flow field affects every other 
part. Many non-linear elliptic equations are solved by 
appropriate relaxation iteration techniques by casting the 
equation in Poisson's form with the non-linearity acting as 
the driving force. Solutions to hyperbolic equations are 


usually obtained by the method of characteristics or by 


finite difference marching techniques which use an artificial 
viscosity to represent the average jump conditions across the 
shock wave. In mixed supersonic and subsonic flows, normal 
numerical procedures break down because the boundary between 
the two regions is not known a priori. 

Finite element numerical techniques have evolved as a 
powerful tool in obtaining approximate solutions to a wide 
variety of engineering problems, particularly ones with 
Neumann-Dirichlet boundary conditions, i.e., elliptical prob- 
lems. They offer several outstanding advantages. Some of 
these are: 


1) Non-homogeneous problems may be treated with relative 
simplicity. 


2) Complex geometries may be modeled with relative ease 
since the elements can be graded in size and shape to 
follow boundaries of arbitrary shape. 

3) Once the finite element model has been established, 

a variety of problems can be solved by supplying the 
computer with the appropriate data. 

Chan and Brashears [Ref. 5] developed a finite element 
computer program for steady transonic flow over a non-lifting 
airfoil. This program uses the least squares method of 
weighted residuals to approximate equation 1 by a system of 
algebraic equations, and assembles the equations in a special 
way to account for the hyperbolic region of flow. This tech- 
nique prevents the influence of downwind nodes from propagating 
upstream in the supersonic region. 

The purpose of this thesis is to investigate the possi- 


bility. of speeding the convergence of a solution by transform- 


ing the steady transonic equation to the unsteady equation 


te 3 
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and integrating over time until the time dependent terms 


vanish and to extend the program of Ref. [5] to the tran- 


sonic region of a converging-diverging nozzle. 


it. DISCUSSION OF THE FINITE ELEMENT APPROACH 


In a continuum problem of any dimension, the field 
variable, whether it is velocity potential, velocity, temper- 
ature, displacement or some other quantity, possesses infin- 
itely many values because it is a function of each generic 
point in the solution region. Consequently, the problem is 
one of an infinite number of unknowns. The finite element 


approach subdivides the solution domain into a finite number 


of subdomains called elements and expresses the unknown field 
variable in terms of assumed approximating functions within 
each element. The approximating functions are sometimes 
called interpolating functions and are defined in terms of 
the values of the field variables at specified points called 
nodes or nodal points. Nodes usually lie on the element 
boundaries where adjacent elements are considered to be con- 
nected. In addition to boundary nodes, an element may also 
contain interior nodes. The nodal values of the field variable 
and the interpolating function for the elements completely . 
define the bahavior of the field variable within the elements. 
Once these unknowns are found, the interpolating functions 
define the field variable throughout the assemblage of the 
element. 

Clearly, the nature of the solution and the degree of 
accuracy of the approximation depends not only on the number 


and size of the elements used but also on the interpolating 
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functions selected. Interpolating functions may not be 
chosen arbitrarily because certain compatibility conditions 
must be satisfied. Often such functions are selected so 
that the field variable or its derivatives are continuous 
across adjoining element boundaries. Once the problem is 
formulated in terms of individual elements, the contributions 
of each element may be assembled to define the entire solu- 
tion domain. This means, for example, that if we are treating 
a problem in stress analysis we can find the force-displace- 
ment or stiffness characteristics of each element and then 
assemble the elements to determine the stiffness of the whole 
structure. Finite element solutions are not, of course, re- 
stricted to structures problems, but the matrix of equations 
defined by the interpolating functions and the nodal field 
variables is still referred to as the stiffness matrix 
regardless of the field variable in the problem. 

Solutions to continuous problems by the finite element 
approach always follow a systematic step-by-step process. 
This process is completely general to the finite element 


method and it is outlined below. [Ref. 4] 


1. Discretize the continuum. 


The first step is to divide the solution domain into 
elements. A variety of element shapes may be used and one 
or more different element shapes may be used in the same 
region. The type and number of the elements used in a given 


problem are a matter of engineering judgement. 
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2. Select the interpolating functions. 


The next step is to choose the type of interpolating 
function to represent the variation of the field variable 
over the element. The field variable may be a scalar, a 
vector, or a higher order tensor. Often polynomials are 
selected as interpolating functions for they are easy to 
integrate and differentiate. The degree of polynomial 
chosen depends on the number of nodes and the nature and 
number of the unknowns and the continuity requirements 
imposed at the nodes and the element boundaries. The un- 
known quantities at the nodes may be assigned to the field 


variable and their derivatives. 


3. Find che element properties. 


After the elements and their interpolating functions 
have been selected, the matrix of algebraic equations which 
express the properties of the individual elements must be 
determined. Several methods are available for this task. 
These are: the direct approach, the variational approach, 
the method of weighted residuals and the energy balance 
approach. Reference [4] is a good source of information 


on the various techniques. 


4. Assemble the element properties to obtain the system 


equations. 


To find the properties of the over-all system, the matrix 


equations expressing the behavior of each element must be 
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added to the matrix equation of all other elements. The 
basis for this assembly procedure stems from the fact that 
connecting elements have common nodes and the field variable 
must be the same for each element sharing that node. 

At this point the boundary conditions for the system of 
equations must be accounted for and the system of equations 


must be modified before it is ready for solution. 


5. Solve the system of equations. 


The assembly process of step 4 produces a set of simul- 
taneous equations which can be solved to obtain the unknown 
field variables. Linear equations have a number of standard 
solution techniques readily available, but solutions to non- 


linear equations are more difficult to obtain. 


6. Make additional computations if desired. 


Important parameters, such as pressure coefficient in 
aerodynamic problems, may now be calculated from the values 


of the field variables. 


a 
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III. THEORY AND BASIC EQUATIONS 


A. STEADY TRANSONIC FLOW 

Chan et al. [Ref. 5] developed an algorithm to analyze 
steady transonic flow over non-lifting thin airfoils. Bound- 
ary layer effects were ignored and the imbedded shock wave 
was assumed to be weak. These assumptions are consistent 
with transonic small disturbance theory which can be expressed 
mathematically by the following expressions. 
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(1 - MZ - ME(y + 1)6,)0,4 + byy = 0 (1) 


Boundary conditions - 


Veo = 0 at infinity (2) 

v = (1 + ujdg/dx on the airfoil (3) 

u=0- on the line of symmetry 
(4) 


where g(x,y) defines the airfoil and dg/dx describes the 
airfoil slope. 

The above expressions appear in their dimensionless form 
where > = perturbed velocity potential and the perturbed é 
velocity components in the x and y directions are respectively 


defined as 


es 
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M, = freestream Mach number and y = ratio of specific heats 
which for air is taken to be 1.4. The physical coordinates 
x' and y' and the velocity potential ¢' are related to the 
dimensionless quantities by 
x « 2 fe, y = "fe, o = o'/cUL 

where c is the chordlength of the airfoil and U. is the free- 
stream velocity. 

Once the flowfield solution is determined in terms of the 


perturbed velocities, the secondary unknowns are subsequently 


calculated. These include: 


U 
a molten <5) 4) gy*t 7 (5) 
V 
= M= a (6) 
gh i aa ital nia aa, ; 
%o =a ts 
2 2 
q+ -[#+ a 3-5] 2 - (8) 


In the above, U, = 1 is the normalized freestream velocity, 
a = local sound speed, p = local static pressure, M = local 
Mach number, V = total velocity, Py = stagnation pressure 


and Cy = pressure coefficient. 


B. UNSTEADY TRANSONIC FLOW 
Unsteady transonic inviscid flow may be expressed in terms 
of the velocity potential ¢(x,y,t) to a first order approxima- 


tion by 
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(1-M2-M2(y +1) 4) 0, +o, - Mo, - M70, = 0 (9) 


yy 
which has the same non-linear coefficient retained for steady 
transonic flow in Equation l. 

Boundary conditions require that the disturbances vanish 


far from the airfoil, 


o, = 0 6% 0 at infinity (10) 


and that the flow remain attached to the body. Let B(x,y,t) =0 
define the body at any instant. The surface tangency restraint 
may now be expressed by the substantial derivitive DB/DT 


vanishing. 


ps/DT = 5, + (1 +9, )B, + OLB. (11) 


If the body remains stationary By = 0, and the tangency 


condition becomes the same as in steady flow, namely 
v = (1 + u)dg/dx (12) 
where dg/dx represents the airfoil slope. 


The pressure coefficient for isentropic unsteady com- 


pressible flow is defined by 


iim «f 3 2 2,,¥/(¥-1 
Cp ae {GM ge tae ye gl a 


Expanding by the binomial expansion and retaining only the 


i 
t 
first order terms gives, 


C, = -20, - 26 (14) 


p t 
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IV. FINITE ELEMENT FORMULATION 


The method of weighted residuals is a technique for 
approximating solutions to linear or non-linear partial dif- 
ferential equations and it is the basis for the finite element 
formulation of the transonic small disturbance equation 
(Equation 1). This procedure involves assuming the general 
functional behavior of the field variable which would approx- 
imately satisfy the basic equation and boundary conditions. 
Substituting this approximation into the original differential 
equation results in some error called a residual. A system 
of algebraic equations results when a weighted average of the 
residual is forced to vanish as it is averaged over the entire 


domain. 


A. STEADY FLOW 
The approximate solution to equation 1 is assumed to be 
a Nid; (15) 
: ; where N; are the interpolating functions which exhibits the 


behavior of equation 1 and $; are the undetermined parameters 


at the nodal points. 


ew: 


When : is substituted into equation 1, the resulting 


residual is 


| R= (1 - M2 MZ + LING INS xe NS yy (26) 
7) 
| } 18 


rem = 


The weighted average is determined by multiplying the 
residual R by m linearly independent weighting functions W; 
and integrating over the elemental domain. Forcing this 


residual to vanish yields, 
SW; RdA = 0 i=iltom 


Chan et al. [Ref. 5] found that when the weighting function 
W; for equation 1 is chosen to be dR/ 39; the resulting matrix 
is positive definite and well conditioned. This choice of 
weighting functions is referred to as the method of least 
squares because it is equivalent to minimizing the square of 
the residuals summed over the domain with respect to the un- 


determined parameters. That is, 
X = rrR2da 
ax/a % = ssaR/ao,RdA = 0 (17) 


Integrating over the domain produces the system of 
algebraic equations 


ick. = 8 
$55; 0 (18) 
where the elemental matrix Si; is defined as 


Sij 7: SIQ5P 5 dA (19) 


With Q; and Ps equal to 


= - 2. - 
Qj 7 C2 - Mo MaGy * DINy OK) NG scx * NG yy 


i 2 
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B. UNSTEADY FLOW 


Development of the unsteady flow finite element equations 
is similar to the procedure used to formulate the finite ele- 
ment equations for steady transonic flow. The least squares 
method of weighted residuals is again used but » is now a 
function of time as well as the spatial coordinates x and y. 


The approximate solution has the form, 
@ = Ni (x,y), (t) (20) 


Substituting $ in the unsteady transonic small disturbance 


equation, the weighted residual becomes, 


x = S(Ry + Ry + Rg) 7dA (21) 
where 
oe Tih IG le IN. EN. t8: 
1 0 30 k k,x°  j,Xx aX" * 3 
= - 2 4 
Ry = -2MN, 6, 


Rs = -MN,S, 
Expanding equation 21 yields 

x = ssUR + R2 + RP + 2R,R, + 2R,Rz]dA (22) 
and minimizing with respect to the undetermined parameters 
5 the following system of algebraic equations result, 

5 * 0 = 2fF8R5/ 995 (Ry + Ry + RgJdA 


aR / 99, sy Ps 


20 
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where P; has been previously defined in the steady finite 


element formulation 


The above equation may be rewritten in the form 
Sty t Se + OA 
ij*j ij*j ij") As} 


The stiffness matrix Si is the same as that developed 


for the steady transonic equation and the damping (SC; 5) and 


mass (SM, 5) matrices are defined below. 


2 | 
SC; = -SIMIN, ,P ida (24) 
meee 
SM; 5 = -SIMIN;P;da (25) 
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V. ELEMENT DESCRIPTION AND ASSEMBLY OF EQUATIONS 


A. ELEMENT DESCRIPTION 

The basic element used in the finite element program is 
the non-conforming cubic triangular element developed by 
Bazeley et al. [Ref. 2]. Also used in the program are the 
quadrilateral and trapezoidal elements constructed from these 
triangular elements. These three types of elements can be 
mixed and used freely in the entire flow region except that 
only trapezoids should be used to cover the supersonic and 
mixed region in order to enact the special assembly procedures 
required by the hyperbolic equation which describes the flow 
in that region. 

The basic triangular element is shown in Fig. 1, which at 
each vertex has the velocity potential and the velocity com- 
ponents as the undetermined parameters. This type of element 
was chosen because both Dirichlet and Neumann boundary condi- 
tions can be imposed with equal convenience. In addition, 
because velocity components are used as primary unknowns 
secondary parameters, such as Mach number and pressure coef- 
ficient can be calculated directly without resorting to 
numerical differentiation, which would produce additional 
errors. 


In the element, the approximate solution is assumed as 


¢ = N.6, (i = 1 to 9) 


22 
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1 (Pi 4.) 


Figure 1 - Triangular Element 
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In which >;'s are the nine undetermined parameters of » and 
N; are the interpolation functions which are expressed in 
terms of the area coordinates. 

The shape or interpolating functions and their first and 
second derivitives are defined below. 


Letting, 


A = area of triangle 1-2-3 = (bj cy - bye 5)/2 


a = 0.5 (cy - ¢,) 
B= 0.5 (bj - by) _ 
H 


= 54555, 
Hy = Diozoy t Djoyeg + OLS 5S; 
Hy = 5S 5S F CpSqS i * ORES; 
Bow * 2(z bj) + 65 >,by + 6b id5) 
Hy = 2(c 55S, + 5 5oKo4 + KC 4o5) 
with i = (1,2,3),% = (3,1,2) then one has for 1 = (1,4,7), 
i= (1,2,3) 
Ny > 67 (3 - 2¢,) + 2H 
Ny x * (6b;c; (1 - 63) + 2H] /20 
N 


1,y * (Sey, G - 64) * 2HJ/24 


24 
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=n ore 
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es [6b; (1- 2¢,) + 2H] /(2a)? 


[6c2 


, (> 2%,) + 2H] / (2h) 2 


Livy 
for 1 = (2,5,8), 1 = (1,2,3) 


= re ' 


om 4 Z 
Ny x [2b &5 (45, C; 4) * 255 + OH] /24 
Nyy 5 [2c, 55 (q 5, - C5 Sy) : oH, | /2 4 
FS 2 7 2 
Na xx = [ 2b; (455 C5 Sy) + 4b; (24)5, + on] /(24) 


a 2 p 2 
Ni yy [2c; (q 5, ©5534) + OH ay] /(24) 
foe 2 «= (356,939), k= 122,3) 


Zz 


a 
f 


= 
u 


[2b; 55 (b; o) - by 5) + i /24 


= - 2 
N (2c; 5; (b; Sy bo 5) + 255 + BH, ] /24 


Z 2 : 2 
Ny xx * (2b; (b; Sy by 55) ¢ BH] / (24) 

“ 2 z 2 
My wy [2c; (b5 54 by) + 4c,(24); + BH] ptony*, 
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Quadrilateral and trapezoidal elements as shown in Fig. 2 | 
are also used in the present program, the former is used in 
the subsonic region and the latter in the mixed and super- | 
sonic region. The element matrix for the quadrilateral ele- . : 
ment is obtained by combining appropriately the matrices for 
two triangles, while the matrix for trapezoidal element is 
obtained by averaging contributions from two left-running and 
two right-running triangles. The averaging process is designed 


to remove the bias effects inherent in the quadrilaterals used. 


B. ASSEMBLY OF EQUATIONS 

Straightforward application of the finite element assembly 
technique to transonic flows would fail (the solution diverges) 
because this would allow disturbances to propagate upwind in 
the supersonic region of flow where the governing equation is 
hyperbolic. Hyperbolic equations have a time-like dependency 
in that the solution at the downwind station is affected by 
the upwind station but not vice-versa. Assembly techniques 
for a transonic flow finite element program must take into 
account this time-like dependency. If the x-axis is taken 
as a time-like direction in the supersonic region, the ele- 
ment matrix may be assembled in a way similar to a backward 
finite difference operator, which has been successful in 
solving hyperbolic equations. 

Consider the rectangular element sketched below with the 
upwind station I and the downwind station II, each having two 
nodal points with the element type chosen. The element 


matrices can be constructed in the usual manner. ‘ 
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a. Quadrilateral 
Element 


b. Trapezoidal 
Element 


Figure 2 - Quadrilateral and Trapezoidal Elements 
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' 
However, before assembling the element matrix into the system 
matrix the non-linear coefficient of equation 1 is evaluated. 


c=1-M-M? wy +i1)u 


The sign of the coefficient being positive, zero, or 
negative defines the equation as elliptic, parabolic, or 
hyperbolic. If C is non-positive for all nodes in the element, 
the rows representing the improper downwind influence on the 
solution at an upwind station are ignored during assembly. 
This feature is automatically applied in the program requiring 
only a little care in arrangement of the nodes of the element. 
In the anticipated supersonic region, element node points 
should be arranged in the order as indicated in figure 3, 
starting with the upper left corner and proceeding in the 
counter-clockwise direction. In the elliptic region, i.e., 
where the coefficient is positive, no special assembly tech- 


nique is invoked. 


C. ITERATIVE PROCEDURES 
With the equations assembled and the proper boundary con- 
ditions imposed, the system of non-linear algebraic equations 
is solved by iterative procedures in the form 
a AR} 6 


th 


to solve for the solution in the n iteration. The function 


o is defined as 


- (n-1) 


o = Oo (n-1) 


+ (1-6) 
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Figure 3 - Nodal Arrangement for Supersonic Region 
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in which the under-relaxation factor 9 is in the range 


‘ 


Q0<68<1. For subsonic flow § = 1, which is simply a succes- 
sive approximation, yields good results, but it is necessary 
to under-relax somewhat with 8 approximately .5 for super- 
critical flow. Generally, a smaller relaxation factor will 
make the solution more stable but it will tend to slow down 
the rate of convergence. 

Equation 23 is subject to the convergence criterion that 
the change in local Mach number between two successive iter- 
ations is less than a prescribed value ¢€ at all nodes in the 
flow field. That is, 

7) 2 y(n) 
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VI. INTEGRATION OF UNSTEADY FINITE ELEMENT EQUATIONS 


The unsteady transonic small disturbance equation (equation 
9) when suitably reduced to a finite element approximation 


appears in the forn, 


One e me te + ee ; 26 
$559; SC; 59; sia = de R; (26) 


This equation is analogous to a damped spring mass system, 


hence §S SC.., and SMi5 are respectively referred to as 


ia” ij 
the stiffness, damping, and mass matrices. 

Mathematically, equation 26 represents a system of second 
order differential equations with constant coefficients, which 
can be solved by standard numerical procedures for differential 
equations, such as Runge-Kutta or Milne methods. However, 
this would be a very costly technique if the coefficient 
Matrices are very large. In practical finite element analysis 
there are a few effective methods which take advantage of the 
banded matrices usually encountered in finite elements. One 
such method is the direct numerical integration method. 

Direct integration involves a numerical step-by-step 
procedure aimed at satisfying equation 26 only at discrete 
time intervals At apart and not over all time t. Conceptually, 
direct integration is a finite element method in space and a 
finite difference method in time. Examples of direct inte- 
gration are the central difference method, Houbolt integration 


and the Wilson method. The first two schemes are finite 


ee 


difference schemes whereas the latter is a linear accelera- 
tion method. Linear acceleration integration assumes a 
linear variation of acceleration from time t to time t+At. 

Central differencing can be very effective in the solu- 
tion of many dynamic problems especially those that involve 
a large system of equations. However, this method is unstable 
for all time steps larger than a critical time step. 

Houbolt integration is an implicit finite differencing 
method related to central differencing, only it has the ad- 
vantage of being stable for all TIME STEPS. 

The Houbolt method was used to integrate the unsteady 
finite element equations because of this stability. 

Houbolt integration uses the following finite difference 


expansions: 


: ; 2 2 
bistent ~ (2Oi tage” Oat * MPile-at” Si, e-2aey)/4t 


Oi steat ~ Cleesae - 18O5 4 * OG bene 2Fi te ane] /Odt 


which are two backward-difference formulas with errors of 


order (at) ?. 


The solution of 9. must satisfy equation 26 and at 


2b btAt 
time tt+tAt equation 26 becomes 


Sipej,teac® Sipe tone * Mijej teat 7 


Substituting the finite difference formulas for F and 


»tr+dt 
rearranging all known vectors on the right hand side, the 


solution for $. 


j t+atis obtained, namely: 
’ . 
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(a,SM; ; + a 


) 
j 19°; * 84) 


j * Sipe jy teat ~ Ryton (27) 


+ (a,SM; + azSC 


j + (a,5M, ; + a-SC. 


ij)?j,t 5905 5)5 t-2at 


Where the constant integration coefficients are: 
a. = 2/at? 

a, = 11/6,at 

a, = S/at 

az = 3/at 

a, = 2a5 

ac = -a,/2 

ap = ay/2 

a, = a,/9 
Equation 27 may be written as 


SE..¢ RE. (28) 


Sd es ae 
where the effective stiffness matrix SE; ; and the effective 


load vector RE, are defined as: 


SEs; = $i; + a)5M;; + a,5C;, 


mg eg artsk” Matyi secge * “0*y, 2-290? 


* SCs Cage ce * Mghy teat * Ore c- zee? 


| Accurate knowledge of the vectors °F teat and *5 6 t-2at 
are required to yield an accurate solution for 2S teat and 
4 
: 3 
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normally the Houbolt integration scheme requires a special 
starting procedure to determine the initial two vectors 

Oj ,at 29d >; one: 
this problem is to integrate the equations until they con- 


However, since the primary interest of 


verge to a steady state solution, it is not necessary to 
obtain an accurate time history of the flow. Errors induced 
by the inaccurate starting vectors will vanish as time 
approaches infinity. Therefore, the starting vectors may be 


chosen somewhat arbitrarily. 
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VII. CONVERGING-DIVERGING NOZZLE 


S. F. Shen [Ref. 10] demonstrated the feasibility of 
calculating compressible flows through a converging-diverging 
(Laval) nozzle by dividing the region of calculations into 
three patches, a subsonic region, a supersonic region and a 
transonic one, of course bounded by the other two regions. 
The locations of the boundaries for each region were chosen 
arbitrarily provided the sonic line is bracketed by the 
subsonic and supersonic boundaries. 

Two different finite element formulations were used for 
the subsonic and supersonic regions, but Shen [10] resorted 
to analytical approximations to cover the transonic patch. 
This restricted the calculations to nozzles with small throat 


curvatures because no analytic solutions exist for nozzles 


with large throat curvatures. It is conceivable that STRANL-II 


could be adapted to provide a continuous solution throughout 
all three regions. 

Outside, the transonic region of flow the governing small 
disturbance equation is 


Z 
(1 - MO, * by * 0 (29) 


This holds for both subsonic and supersonic flow. Comparing 
equation 29 with equation 1, the transonic small disturbance 
equation, we notice that only the non-linear coefficient 


M2 Cy + l)u distinguishes the two equations from each other. 
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This coefficient becomes negligible when the Mach number 
becomes less than .8 or greater than 1.2. With this con- 
sideration in mind, it was assumed that equation 1 would 
adequately describe the flow through the nozzle and that 
the finite element formulation developed for the non-lifting 


airfoil would apply to the Laval nozzle. 


A. BOUNDARY CONDITIONS 

Two solutions are possible for a converging-diverging 
nozzle: 1) Symmetric flow, where the flow is subsonic through 
the domain, except for a small supersonic region near the wall 
in the throat, and 2) Asymmetric solution, where the flow 
accelerates to sonic velocity in the throat and then continues 
to accelerate to supersonic velocity in the diverging section. 
Different boundary conditions apply for the two solutions. 
For the symmetric case, both inlet and exit velocities must 
be specified. Inlet and exit velocities are equivalent in 
the subsonic solution. The supersonic solution requires that 
only the inlet velocities be specified. If the exit veloci- 
ties are also applied, the problem is overspecified and the 
solution may not converge. 

Velocities at the inlet and exit are not uniform in the 
y direction, therefore the disturbances cannot be set to zero 
as in the case of the non-lifting airfoil. Boundary velocities 
must be calculated by solving equation 29 analytically. 

Equation 29 is a linear equation which can be mapped to 
Laplace's equation, 


vos 0, 
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by letting y' =,/ (1 - M,) y. Laplace's equation is easily 


solved for the case of the hyperbolic nozzle by transforming 
from cartesian coordinates to elliptic coordinates. This 
transformation simplifies the solution because the stream 
lines must be hyperbolas to follow the nozzle boundary and 
therefore follow the hyperbolic coordinate y = constant. 

If the elliptic coordinates , and y are chosen such that 
the curves » = constant are ellipses and the v curves are 
hyperbolas, then the velocity potential which satisfies 


Laplace's equation for a hyperbolic nozzle is simply 
where A.is a constant of integration. The stream function is 


The transformation w = u + iv = cosh 4(2z/a) gives rise 


to the elliptic coordinates 


y = 1/2 a cos'.ucos v, x = 1/2 a sinhusin v 


£ L/2 @ (cosh - sin*v] 


t;.:* cy + a/2)° + x? 


tT, * (y - a/2)? + x? 


Solving for uy and v produces 


1 


u = cosh” [(r, + r,)/a] 


1 


v = cos” ((r, + r,)/a] 


The nozzle boundary is defined by Vy = constant, which 


along with the equation for the nozzle wall in cartesian 


coordinates, y' = f(x), implicitly defines the constant a. 


Substituting for »p in the velocity potential produces, 


A cosh”+ ((r, + r,)/a] 


from which the velocities may be determined. 


= = + — 
Pe Ait, * Ealfay> = 1 Gy Sy 


; or or 
a ee ee Shakey regen 2 
¥ ee ee dy dy 
ct+——*)* 1 
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x 
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fly + a/3)% + x2 
9,2 x 

ox afly * a/>)¢ + x? 
ory Y #-a/Z 

dy Pr 5 #2)” ae 
or2 Piers 


ay Ky - a/,)2 + x2 


The constant of iutegration A may be determined by speci- 
fying the flow rate through the nozzle, but when the inlet 
velocities are normalized with respect to the freestream 
velocity (U,) A is factored out of the problem. 

Velocities for compressible flow can be solved by mapping 
back to the physical coordinate system (x,y plane). 

Other boundary conditions are universal to both problems. 
These are: 

u = 0 on the line of symmetry 


v = (1+ u)df/dx at the nozzle wall 


A 


F(x) defines the nozzle boundary in terms of a ratio of the 
throat semi-height as a function of x. The throat semi-height 
is taken to be 1 for convenience. . 

Pressure ratio, sound speed, and Mach number are calculated 


’ ‘ 


as before by equations 5 through at 


, 
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VIII. DISCUSSION OF RESULTS 


A. TIME INTEGRATION TO A STEADY STATE SOLUTION 

As stated before, the Houbolt method of integration is 
stable for all time steps. Results of the test cases bear 
this out with the larger time steps providing the most rapid 
convergence to a steady state solution. Time steps were 
tried from t = .1 to t = 100. Time steps larger than t = 100 
were not attempted because as t becomes too large the influ- 
ence that the damping and mass matrices have on the effective 
stiffness matrix becomes negligible compared to the stiffness 


Matrix. That is* 


Zz 
SEs; Si; 2 28C 5 ;/at + 6SM, ,/At 


as At + @ 


SE i; ~ Sij 


The starting solutions were chosen somewhat arbitrarily. 
$ ; zat Was chosen to satisfy the first iteration of the steady 
~? 
solution 
Sia cds = 
ij?j,2at 
when the non-linear term (u) in the coefficient 


1- M2 - M@(y+1)u 


was set to zero. $:; and 9. were chosen as multiples of 
jJ,dt j,0 
% 5, 2at and respectively they were 


This starting procedure proved to be superior to chosing 
the first three vectors closer to the converged solution. 


If 9 and 9 were chosen to be the last three 


j,2dt? °j,ar? j,0 

time steps of the previous case, the solution oscillated and 
converged much slower than with the starting solutions chosen 
as above. 

The stiffness, mass, and damping matrices were recalculated 
after each time step, using the under-relaxation technique 
described above. This was necessary to utilize the special 
assembly procedures invoked by STRANL-II to prevent the in- 
admissable influence of downwind nodes from propagating up- 
stream in the supersonic region. 

For barely critical flow (M,= .861) and subsonic flow, 
an under-relaxation factor 9 = 1 (successive approximation) 


resulted in convergence to a steady state solution after only 


three time steps. Eleven time steps were required for the 


supercritical solution to converge using the same relaxation is 
factor. Reducing § to .5 increased the rate of convergence 

and the solution achieved steady state after six time steps. 
Figure 4 compares the steady state solution for a 6% thick 
circular arc airfoil at M, = .909, using the same integration 
method, with the results obtained in Ref. [5]. Chan's results 
converged in 10 iterations after using the results from the 


barely critical flow as an initial guess to the 
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supercritical solution. Figure 4 is a plot of local Mach 


numbers at boundary nodes on the airfoil. 


B. CONVERGING-DIVERGING NOZZLE 
The nozzle chosen for the test cases was the two-dimen- 


sional Oswatitsch nozzle with the boundary defined as 


Fa Oe ie 28 


where the throat at x = 2.5 has a semi-height of 1. The inlet 
was taken to be x = 0 and the exit was at x = 5. M, = .44, 
the inlet Mach number on the nozzle center-line was chosen to 
yield sonic conditions in the throat. 

Two solutions were possible for this inlet condition,- the 
symmetric solution and the asymmetric solution; but neither 
solution was achieved by the finite element method. Although 
the solution converged for the subsonic case in three itera- 
tions, center-line Mach numbers deviated significantly from 
both one-dimensional theory and from Oswatitsch's approxima- 
tion (Fig. 5]. When the local Mach number M exceeded the 
inlet Mach number by approximately .2 (M > .64) the solution 
was invalid. Differences at the center part of the nozzle 
are due to an essentially incorrect free stream Mach number. 
Patching the solution at x ~*~ 1.5 would improve the solution. 

A second test case was run for the supersonic section of 
the nozzle with the inlet boundary on the sonic line. The 


exit boundary was left free to float. Here the solution was 


unstable and no meaningful results were obtained. 
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IX. PROGRAM MODIFICATION 


The finite element computer program for non-lifting air- 
voils, as developed in Ref. [5], is separated into two parts. 
These have been designated STRANL-I and STRANL-II by Lockheed 
Corporation. STRANL-I generates a finite element mesh to be 
used as inputs to STRANL-II, which assembles the finite ele- 
ment equations, applies the boundary conditions, and solves 
the non-linear system of equations. Detailed descriptions 
and instructions for the use of the two programs can be found 
in Ref. [5]. Only the modifications to the above programs 


will be discussed in this section. 


A. UNSTEADY EQUATIONS 
Modifications to STRANL-II to form and solve the unsteady 
finite element equations were three-fold: 


1) The new elemental matrices SCj 


and SM jj were calculated 
and assembled. 


j 

2) All the matrices were stored on an external magnetic 
disk to be accessed and reassembled later because of 
the amount of space required to store three large 
matrices, in core memory. 

3) The effective stiffness matrix SE;; and the effective 
load vector RE;; were assembled, and the system of 
equations solved. 

Several existing subroutines in the original STRANL-II 
program were modified to assemble the damping and mass matrices. 
These include subroutines NEWK, EMTC, DERV, and EMQT. Two new 


subroutines were added to perform the other tasks. 


4$ 


TOD 


B. MODIFICATIONS TO CALCULATE THE MASS AND DAMPING MATRICES 
EMTC in the STRANL-II program calculated the elemental 


stiffness matrix by numerically integrating the equation, 


| Si; =  Q,P;da 


Equations were added to EMTC to perform the additional numer- 
ical integrations for the mass and damping matrices. All 
three matrices Were calculated at the same time. EMQT 
assembled the elemental stiffness matrices for a quadrilateral 
and trapezoidal element from the contributions of the tri- 
angular elements. Mass and damping matrices were calculated 
in the same fashion. 
| Subroutine NEWK, an original subroutine in STRANL-II, 
which assembled the finite. element equations for steady flow, 
was modified to assemble SC; and SM 5: A new calling argu- 
ment, NMAT was passed to NEWK, which assembled contributions 
from the triangular, quadrilateral and trapezoidal element 


SC... 


Matrices into the global matrices S; ij? 


and SM; 


a? ad 


depending on NMAT being 1, 2 or 3. 
1. Subroutine STORE 


/ 
Given a non-symmetric matrix stored in a banded node, 
plus the right hand side vector, subroutine STORE separates 


this system into two matrices and stores them on a magnetic 


disk. Figures 6 and 7 show the decomposition of a banded 
matrix into banded storage, and the further decomposition 
of this banded stored matrix to two smaller matrices by sub- 


routine STORE. In these figures, D, L, and U represent the 
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diagonal matrix, the lower triangular matrix, and the upper 
triangular matrix respectively. HBW is the half bandwidth 
and R is the right hand side vector. 

STORE requires an additional work area one-half the 
size of the originally dimensioned matrix which is to be 
stored. 

2. Subroutine TIME 


Subroutine TIME integrates the system 


*. * 6¢,.4, + Sas, = kh. 
eae Dae 5 mia” j 


by Houbolt integration. 

TIME reassembles the three matrices which were stored 
on the magnetic disk to form the effective stiffness matrix 
and the effective load vector. Once this system of equations 
is assembled, a banded equation solver is called to yield the 
solution for 9. 


j,t+At’ 
TIME. In Fig. 8 when L=1, the lower triangular matrix and 


Figure 8 is a schematic flow chart of 


the diagonal of the effective stiffness matrix are formed by 
adding the appropriate contributions from the stiffness, mass 
and damping matrices. When L=2, the upper triangular matrix 


is formed in like fashion. 


C. CONVERGING-DIVERGING NOZZLE 
1. Application of the Boundary Conditions 
Regardless of the type of problem for which a set of 
system equations have been assembled, the equations will have 


the form 


READ WORK 


+ UM X WORK 


SOLVE FOR 


Peeat 


RETURN 


Figure 8 - Flow Chart of Subroutine TIME 
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Figure 9 - Flow Chart of STANL-II Modification 
to Integrate Unsteady Equations 
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in which Ki; is annxn matrix and xX; and R; are vectors of 


length n. These equations do not take into account the known 
values of X; on the boundaries. However, for a unique solu- 
tion of the above equation, at least one or more nodal var- 
iables must be specified and Ki; must be modified to render 
it non-singular. For each equation i, either X; or R; must 
be specified but it is physically impossible to specify both 
Xi and R;. There are a number of ways to apply the boundary 
conditions to the equations and when they are applied the 
number of equations is reduced. However, it is convenient 

to leave the number of equations unchanged to avoid major 


restructuring of the computer storage. One such method is 


described below. 


é 


If k is the subscript of the prescribed nodal variable, 


th th 


the k row and the k column of the original K.. matrix are 


set to zero, Kix is set to 1 and Ry 1s replaced * the known 
value of X;,. Each of the n-1l remaining terms of R; is modi- 
fied by subtracting from it the value KiyXpe This procedure 
is repeated for all the boundary values. Of course, when the 
matrix is stgred in a banded mode, the algorithm will differ 
from that for the nxn square matrix, but the procedure is 
similar. 

Subroutine BNDRY applies the boundary conditions for 
the modified program. Setting the option parameter IOPT(4), 
in STRANL-II, equal to 1 will call BNDRY which will read the 


boundary velocities and apply them to a banded stored matrix 


by the method described above. 
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| When the value of the xX, is zero, as in the non- 
lifting airfoil problem, the algorithm becomes simpler than 
the above method because there is no need to either set the 


th 


k column to zero or subtract the value Ki %& from the 


right hand vector. 


4 
X. TEST CASES 


In computing the flow field for either the non-lifting 
airfoil or for the Laval nozzle, the following procedures 
were followed: 


1) The desired mesh was sketched with each node assigned 
a number. 


2) Appropriate input cards based on the sketch were pre- 
pared and supplied to STRANL-I to generate the data 
on punched cards for input to STRANL-II. 

3) The above punched cards were supplied to STRANL-II 
with three additional cards as input parameters for 
each case, plus additional cards for the boundary 
velocities, if the nozzle solution is desired. 

4) Results of the finite element caiculations are printed 
after each iteration, and the converged solution is 
punched on cards for possible later use. 

A. TIME INTEGRATION TO A STEADY STATE SOLUTION 

Test cases for the integration of the unsteady transonic 
finite element equations were conducted to calculate the 
steady transonic flow over a 6% thick circular arc airfoil. 

These tests were made using the same airfoil, mesh, and free- 

stream Mach numbers as Chan et al. [Ref. 5] published. These 


conditions were chosen to provide a source for comparison of 


the results. 


Freestream Mach numbers used in these calculations were: 
M, = -806 (subcritical) 
M, = -861 (barely Critical) 


M, = -909 (supercritical). 
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Each case was treated individually with $j = 0 used as 
the initial guess for each case, whereas Chan et al. (Ref. 5] 


used zero for the initial guess for 4. * 806 and then used 


the computed results as the initial guess for each subsequent 


case. 
DELT, the value for the time step, is input by 2 parameter 
specified in columns 41-45 of the second card following the 


title card for each case when the unsteady option (IOPT (6) = 1) 


is selected. 


ie STRANL- I Program 
a. Input 


Input cards used to generate the finite element 


mesh are listed on the next page. Cards were arranged in 


accordance with Ref. [5]> in the following order: 


Title card 

Option card 

Element cards 

Card for the total number of nodes 

Node coordinate cards 

Card for the number of boundary nodes 

Card for the boundary nodes at infinity 
Card for the nodes on the line of symmetry 
Card for the nodes on the airfoil 

Cards for the slope of the airfoil. 


Input cards to STRANL-I for these tests are listed 


on the next three pages. 


b. Output 
Output from STRANL-I is in the form of printouts 


and punched cards. Printouts from STRANL-I are listed on the 


following eight pages. 
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STRANL-II Program 
a. Input 


Listed on the next three pages are the input cards 
to the STANL-II program. 
b. Output 
The output from this program is in the form of 
printouts for each iteration and punched cards for the con- 
verged solution. Output from STRANL-II for M, = .909 is listed 


on the following eight pages. 
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B. CONVERGING-DIVERGING NOZZLE 
Two test cases were run for the converging-diverging 
nozzle. The first case was for symmetric flow designed to 
yield sonic conditions in the throat of the nozzle. The 
second case dealt with supersonic flow in the diverging sec- 
tion by starting with the sonic line as the boundary of the 
nozzle mesh. Oswatitisch's two-dimensional nozzle [Ref. 9], 
y = 1 %£2{x = 2.8)" with a semi-throat height of 1 at 
x = 2.5 was used in both cases. Boundaries for the subsonic 
nozzle were at x = 0 and x = S. 
1. Symmetric Solution 

This problem was analyzed using the mesh shown in 
Fig. 10, which consists of 126 elements and 152 nodes. In 
the second card (the option card) IOPT(4) = 1 indicates that 
non-zero boundary velocities at the inlet and the exit will 
be read and applied by subroutine BNDRY. This option requires 
that the number of inlet and exit boundary nodes be specified 
in columns 36-40 of the next card. Perturbation velocities 
at the boundary nodes follow on the subsequent four cards. 
Subroutine BNDRY reads u and v respectively for the first 
boundary node and then continues reading u and v for each 
inlet and then each exit node in the order specified on the 
appropriate card. 

2. Supersonic Case--Diverging Section 
The mesh used for this case is sketched in Fig. 11 and 


input cards to STRANL-II follow on the next page. For this 


case the options in effect are IOPT(4) = 1 and IOPT(S) = 
which cause non-zero boundary velocities to be applied to 


the sonic line only. 
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Figure 10 - Finite Element Mesh for the Converging-Diverging Nozzle. 
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